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FOURIER PHASE RETRIEVAL WITH A SINGLE MASK BY 
DOUGLAS-RACHFORD ALGORITHM 

PENGWEN CHEN AND ALBERT FANNJIANG 


Abstract. Douglas-Rachford (DR) algorithm is analyzed for Fourier phase retrieval with 
a single random phase mask. Local, geometric convergence to a unique fixed point is proved 
with numerical demonstration of global convergence. 
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1. Introduction 


X-ray crystallography has been the preferred technology for determining the structure of a 
biological molecule over the past hundred years. The method, however, is limited by crystal 
quality, radiation damage and phase determination 46 . The hrst two problems call for 


large crystals that yield sufficient diffraction intensities while reducing the dose to individual 
molecules in the crystal. The difficulty of growing large, well-diffracting crystals is thus the 
major bottleneck of X-ray crystallography - a necessary experimental step that can range 
from merely challenging to pretty much impossible, particularly for large macromolecular 
assemblies and membrane proteins. 

By boosting the brightness of available X-rays by 10 orders of magnitude and producing 
pulses well below 100 fs duration. X-ray free electron lasers (XFEL) offer the possibility of 
extending structural studies to single, non-crystalline particles or molecules by using short 
intense pulses that out-run radiation damage, thus circumventing the hrst two aforemen¬ 


tioned problems 48 . In the so-called diffract-before-destruct approach 19,20,53 , a stream 


of particles is howed across the XFEL beam and randomly hit by a single X-ray pulse, 
forming a single diffraction pattern before being vaporized as a nano-plasma burst. Each 
diffraction pattern contains certain information about the planar projection of the scattering 
potential of the object along the direction of the beam which is to be recovered by phase 
retrieval techniques [lo) . 

The modern approach to phase retrieval for non-periodic objects roughly starts with the 


Gerchberg-Saxton algorithm 32 , followed by its variant. Error Reduction (ER), and the 
more powerful Hybrid-Input-Output (HIO) algorithm pO 
of the standard iterative transform algorithms (ITA) [8,41 


31 . These form the cornerstones 


















Figure 1. Conceptual layout of coherent lensless imaging with a random 
mask (left) before (for random illumination) or (right) behind (for wavefront 
sensing) the object, (middle) The diffraction pattern measured without a mask 
has a larger dynamic range. The color bar is on a logarithmic scale. 


However, the standard ITA tend to stagnate and do not perform well without additional 
prior information, such as tight support and positivity. The reason is that the plain diffrac¬ 
tion pattern alone does not guarantee uniqueness of solution (see 52 , however, for unique¬ 
ness under additional prior information). On the contrary, many phasing solutions exist for 


a given diffraction pattern, resulting in what is called the phase problem 35 


To this end, a promising approach is to measure the diffraction pattern with a single 
random mask and use the coded diffraction pattern as the data. As shown in 27 , the 


uniqueness of solution is restored with a high probability given any scattering potential 
whose value is restricted to a known sector (say, the upper half plane) of the complex plane 


(see Proposition 6.1) 


Indeed, the sector constraint is a practical, realistic condition to impose on almost all 
materials as the imaginary part of the scattering potential is proportional to the (positive) 
extinction coefficient with the upper half plane as the sector constraint 10 . For X-ray 
regime, the real part of the scattering potential is typically slightly negative which with a 


nonnegative imaginary part gives to the second quadrant as the sector constraint 17 


What happens if the sector condition is not met and consequently one coded diffraction 
pattern is not enough to ensure uniqueness? This question is particularly pertinent to the 
diffract-before-destruct approach as the particle can not withstand the radiation damage 
from more than one XFEL pulses. 

A plausible measurement scheme is to guide the transmitted held (the transmission func¬ 
tion [To]) from a planar illumination through a beam splitter ISO], generating two copies 












of the transmitted field which are then measured separately as a coded diffraction pattern 
and a plain diffraction pattern. In this set-up, the object function is the transmitted held 
behind the particle and the phase retrieval problem becomes the wave-front reconstruction 
problem 10,34 . In practice beam splitters and the masks (or any measurement devices) 


should be used as sparingly as possible to avoid introducing excessive measurement noises. 

With two diffraction patterns, the uniqueness of solution (in the above sense) is restored 
almost surely without the sector constraint (see Proposition 6.1 and Remark 6.3). 

With the uniqueness-ensuring sampling schemes (Section 1.1), ad hoc combinations of 


28,29 


members of ITA (such as HIO and ER) can be devised to recover the true solution 
There is, however, no convergence proof for these algorithms. 

The main goal of the paper is to prove the local, geometric convergence of the Douglas- 
Rachford (DR) algorithm to a unique hxed point in the case of one or two oversampled 
diffraction patterns (Theorems 4.1, 5.3 and 6.2) and demonstrate global convergence numer¬ 
ically (Section]^. 

DR has the following general form: Let Pi and P 2 be the projections onto the two con¬ 
straint sets, respectively. For phase retrieval. Pi describes the projection onto the set of 
diffracted helds (instead of diffraction patterns) and P 2 the data htting. The Douglas- 


Rachford (DR) algorithm is dehned by the iteration scheme 24,39 

+Pi (2P2 


(1) + Di(2P2 - /)2/(") - P2y^^\ fc = i, 2, 3... 

Closely related to HIO, DR also belongs to the ITA family (Section]^. ITA are computa¬ 
tionally efficient thanks to the fast Fourier transform (FFT) and explicit nature of Pi,P 2 
(see 0 below). 


1.1. Oversampled diffraction patterns. Next we describe our sampling schemes before 
we can properly introduce Pi,P 2 and the Douglas-Rachford algorithm for phase retrieval 
(Section]^. 

Let /(n) be a discrete object function with n = (ni,n 2 , ■■■,nd) G Z'^. Consider the object 
space consisting of all functions supported in A4 = {0 < mi < Mi,0 < m 2 < M 2 , ...,0 < 
< M^}. We assume d> 2. 

With a coherent illumination under the Fraunhofer approximation, the free-space prop¬ 
agation between the object plane and the sensor plane can be described by the Fourier 
transform (with the proper coordinates and normalization). However, only the intensi¬ 
ties of the Fourier transform are measured on the sensor plane and constitute the so called 
diffraction pattern given by 

M 

/(m + n)7(i^e-*^™-'", w = (w;i,..., w;,) G [0,1]", M = (Mi,..., M,) 

n=—M m+nSAd 

which is the Fourier transform of the autocorrelation 

^/(n) = /(ii^ + n)/(m). 

mSAd 

Here and below the over-line notation means complex conjugacy. 

Note that P/ is dehned on the enlarged grid 

M = {(mi,..., rud) G Z'^ : -Mi < mi < Mi,..., -Md < rud < Md} 
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whose cardinality is roughly 2'^ times that of A4. Hence by sampling the diffraction pattern 
on the grid 

^ f, , , 1 2 2M. ) 

£-^(wi,...,Wd) I wj - 0, + 1 ’ 2M,- + 1 ’ ■■■’ 2M,- + 11 

we can recover the autocorrelation function by the inverse Fourier transform. This is the 
standard oversampling with which the diffraction pattern and the autocorrelation function 
become equivalent via the Fourier transform |4^[44] . 

A coded diffraction pattern is measured with a mask whose effect is multiplicative and 
results in a masked object of the form g{n) = /(n)/i(n) where {/i(n)} is an array of random 
variables representing the mask. In other words, a coded diffraction pattern is just the plain 
diffraction pattern of a masked object. 

We will focus on the effect of random phases 0(n) in the mask function p(n) = |/i|(n)e®'^^“^ 
where 0(n) are independent, continuous real-valued random variables and |/i|(n) ^ 0,Vn G 
M. (i.e. the mask is transparent). 

Without loss of generality, we assume |/i(n)| = 1, Vn which gives rise to a phase mask and 
an isometric propagation matrix 


( 2 ) 


(1-mask ) A* = c$ diag{/r}. 


i.e. AA* = I (with a proper choice of the normalizing constant c), where <F is the oversampled 
d-dimensional discrete Fourier transform (DFT). Specihcally $ G is the sub-column 

matrix of the standard DFT on the extended grid M. where |A1| is the cardinality of AA.. 

When two phase masks p.i,/i 2 are deployed, the propagation matrix A* is the stacked 
coded DFTs, i.e. 


( 3 ) 


(2-mask case) 


A* = c 


$ diag{/ii} 
$ diag{/i 2 } 


With proper normalization. A* is isometric. 

In line with the spirit of simplifying measurement complexity discussed above, we remove 
the second mask (i.e. p 2 = 1) and consider the propagation matrix 

<h diag{/i} 


27-29 


( 4 ) 


(1^-mask case) A* = c 


<F 


normalized to be isometric, where pi is independently and continuously distributed over 
the unit circle. In other words, one oversampled coded pattern and one oversampled plain 
pattern are used for reconstruction. 

For convenience, we shall refer to this set-up as the 14-mask case to distinguish it from 
the one- and two-mask cases. This and its extension to the multi-mask cases would be the 
set-up for all our numerical simulations in Section 


1.2. Related literature. For the optical spectrum, experiments with coded diffraction pat¬ 
terns are not new and can be implemented by computer generated holograms 11 , random 
phase plates and liquid crystal phase-only panels 26 . Recently, a phase mask with 


randomly distributed pinholes has been implemented for soft X-ray 40 


Coded-aperture phase retrieval was formulated as a convex trace-norm minimization prob¬ 
lem in 


12 


15,18 whose uniqueness was proved in 14 under the assumption that the number 


of independently coded diffraction patterns is sufficiently large (polylogarithmic in |W1|). 





















The convex formulation [12|14 ,15 certainly has a tremendous appeal as global convergence 
can be expected for any proper numerical implementations. However, due to the lift to much 
higher dimensions, the convex program may be computationally expensive (see also [7|[^). 

Alternative non-convex minimization formulations were proposed and solved by various 
gradient methods 13,^. In practice, these algorithms are locally convergent with a com¬ 
paratively large number (> 6) of coded diffraction patterns. 

An important difference between the measurement schemes in these papers and the present 
work (as well as |27| - [^ ) is that their coded diffraction patterns are not oversampled. Another 
distinctive feature of the present setting is that the dimension d > 2 is required for the 
spectral gap (Theorem |5. 3 ) and the uniqueness of hxed point (Theorem 6.2), but not for the 
structure theorem (Theorem 4.1). 

In this connection, we emphasize that reducing the number of coded diffraction patterns 
is crucial for the diffract-before-destruct approach and, in our view, oversampling is a small 
price to pay with current sensor technologies. 

A bigger price may be that we lose the robust injectivity property prized in these works 
(also see [^[^). In other words, just one or two random masks, the phase retrieval map F 
defined as F{x) = \A*x\ with A* given by (|^ or ([^ is injective only after a certain finite 
set is excluded from On the other hand, for any given /, with probability one the 

selection of the mask(s) is such that no other object, modulo a phase factor, produces the 
same data as \A*f\ with A* given by (|^ or (|^. Our numerical results show that this notion 
of uniqueness appears suffices for most practical purposes when DR is implemented for phase 
retrieval. 

However, to the best of our knowledge, even local convergence is not known for any ITA for 
Fourier phase retrieval. The present paper aims to fill this gap. While we can not prove global 
convergence of DR as for the convex setting, we will present strong numerical evidence for 
global convergence. In 21 , we prove the local, geometric convergence of the error reduction 

Local convergence 

was recently proved in 


algorithm with the same measurement scheme described in Section |1.1 
of ER in the case of Gaussian matrix A 


47 


There is much more literature on phase retrieval with generic frames and independent 
random matrices [2]-[6, 15,16,22,25,54,56 which have a somewhat different ffavor from that 
with Fourier phase retrieval as the latter is neither of the former two. Consequently, some 
of our results (Theorems 5.3 and 6.2) can not extend to these other cases without major 
revision. The structure theorem (Theorem 4.1) remains valid, however (see Remark 4.4). 

There also is a growing body of work on phase retrieval under sparsity assumptions, 
see and the references therein. 

The rest of the paper is organized as follows. In Section we simplify the notation for 
presenting the main results and the proof of local convergence. In Section we describe 
the DR algorithm widely used in convex optimization problems and formulate it for the 
non-convex problem of phase retrieval. In Section we prove local convergence of FDR 
under the spectral gap assumption (Theorem 4.1). In Section]^ we prove the spectral gap 
condition for any number of oversampled diffraction patterns (Theorem 5.3). In Section]^ 
we prove that the fixed point of DR is unique for one or two oversampled diffraction patterns 
(Theorem 6.2). In Sectionj^we give numerical examples and demonstrate global convergence 
of the DR scheme. 
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2. Set-up and notation 


For the main text below, we shall simplify the notion as follows. The more elaborate 


notion of Section 1.1 will be useful again in the appendix. 


First, we convert the d-dimensional grid into an ordered set of index. The unknown object 
will now be denoted by xq E C”, n = |A1|. In other wor ds, x p is the vectorized version of the 
object function / supported in C Z^, d>2 (Section |l.l[ ). 

Rank-2 property: xq is rank-2 (or higher) if the convex hull of supp{/} is two (or higher) 
dimensional. 


Let df be a nonempty closed convex set in C" and 


(5) 


x\x = arg min llx' — x 

x'&X 


the projection onto X. 


Sector constraint: xq satishes the sector constraint if the principal value of argxo(j), Vj 
is restricted to a sector [— ott,/^tt] C (— 7r,7r],Vn. As mentioned above almost all scattering 
potentials / have a nonnegative imaginary part and hence satisfy the sector constraint with 
a = 0, P = 1. The sector constraint serves as transition between the standard positivity 
constraint {a = (3 = 0) and the null constraint {a = f3 = 1). 

The sector projection is explicitly given as follows: For j < n 


( 6 ) [xUj) 


x{j) 
3f?[x(j)e 


—2/37rl 


0 


if Zx(j) G [—avr,/? 7 r] 
if Zx(j) e [/3ti, {(3 + l/2)7r] 
if Zx(j) G [—(a + l/2)7r, —avr] 
else 


and [x]x{j) = 0 , j > n-f 1 . 


Phase retrieval problem. For a given unknown object xq of rank > 2, let A* = [a*\ G 
0 Axn pg propagation matrix given by 0 . 0 or (|^ where A* is normalized to be 
isometric and b = |A*xo| G be the data vector. Phase retrieval is to find a solution x to 
the equation 

(7) b = |A*x|, X E X. 

We focus on two cases. 

1) One-pattern case: A* is given by ([^, [x];r is given by 0. 

2) Two-pattern case: A* is given by ([^ or (|^, A = C"^. 

Phasing solution is unique only up to a constant of modulus one no matter how many 
coded diffraction patterns are measured. Thus the proper error metric for an estimate x of 
the true solution xq is given by 

( 8 ) min ||e“*®xo — x|| = min ||e*®x — xnll. 

eeiR eeiR 
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Throughout the paper, we assume the canonical embedding 

C*" C C" C C^, n<h<N. 


For example, if a; G C”, then the embedded vector in C” or C^, still denoted by x, has zero 
components x{j) = 0 for j > n + 1. This is referred to as zero padding and h/n is the padding 
ratio. Conversely, if a; G O’* or C^, then [x]„ G denotes the projected vector onto €"■. 
Clearly, [a;]c« = [x]n- 

The vector space C'^ = ©r iR^ is isomorphic to via the map 


(9) 


G{v) : = 


3(n) 


Vn G C 


N 


and endowed with the real inner product 

{u,v) := iR.{u*v) = G{uyG{v), G C'^. 


With a slight abuse of notation, we will use G{u) to denote the conversion of a complex¬ 
valued vector M in €"■, C” or to its real-valued version. 


Phase factor: Let y Q y' and y/y' be the component-wise multiplication and division 
between two vectors r/, y\ respectively. For any y G C'^ dehne the phase vector oj G with 
= l/(j)/ll/(j)l where \y{i)\ 7 ^ 0. When \y{i)\ = 0 the phase can be assigned arbitrarily 
and we set uj{i) = 1 unless otherwise specified. 


3. Douglas-Rachford algorithm 

Phase retrieval can be formulated as the following feasibility problem in the Fourier domain 
(10) Find yeA*Xr}y, y-.= {y e : \y\ = b}. 

Let Pi be the projection onto A*X and P 2 the projection onto 3^: 


( 11 ) 


y 


P,y = A*[AyU, P.,y = b Q ^ 


\y\ 


Then DR (|^ becomes with 


( 12 ) 


Sf{y) = y + A* 


R( 260 1^1 y 


A 


hr^ y 
60 ^ 

\y\ 


which we call the Fourier-domain DR (FDR) to contrast with the following object domain 
version^ 

Let A* = G be a complex isometric extension of A*, implying that A±A]_ = 

/, AA*j_ = Aj_A* = 0. Then the phase retrieval problem can be more generally formulated as 
\A*x\ = b,x & X. Consider the feasibility problem 

(13) Find xexnx, := |x G C” : |i*x| = 6 |. 

Let Pi be the projection onto X, i.e. Pix = [x\x, and P 2 the projection onto X. When 
n = N (hence A is unitary). 


(14) 


P 2 X = A \ bQ 


A*x 

\A*x\ 
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and ( 12 ) is equivalent to 


(16) 


S{x) = X + 


A{2bQ 


A*x 

\A*x\ 


X 


AlbQ 


J A) 


A*x 

\A*x\ 


In this case, we have 

(16) A*SA = Si, for fi = N. 

In the 1-pattern case with the standard oversampling N = n ^ An, A = A is unitary and 


(15) is also known as the Hybrid-Input-Output (HIO) algorithm (with the HIO parameter 
set to one) [8,30 


For h < N (as with two oversampled patterns N ^ 8n with the standard padding n An), 
the precise form of P 2 is not known explicitly. For the purpose of contrasting with (12) and 


for lack of a better term we shall call (15) (with n < N) the generalized Object-domain 


Douglas-Rachford algorithm (ODR for short). The ODR family is an interpolation between 
the HIO and FDR. 

While ODR depends explicitly on h, FDR is independent of h in the sense that 


(17) 


S,(!/)=9 + i- 




A) 


y 

- 60 ^ 

\y\ 


since [Ay]x = [Ay]x G O’* and A*[Ay]x = A*[Ay]x. 


A. Local convergence 


For simplicity, we shall analyze FDR (|12|) without the sector condition: 
(18) 


S,(y) y + A’ A\2b& - y) -i® 


whereas ODR (15) becomes 
(19) S{x) = X + 


H 26 0 


A*x 

\A*x\ 


X 


A 60 


A*x 

\A*x\ 


However, we make no assumption about the number of diffraction patterns which can well 
be one. 

The main result is local, geometric convergence of the FDR algorithm. 


Theorem 4.1. Let Xq G C” and A* any isometric N x n matrix. Suppose N >2n and 

( 20 ) 


max \\u\ 

uGC'^ ' 


I ^||3(i?*M)|| <1, B ■.= A diag * ^ 


\A*xo\ 


Let be an FDR sequence and x^^^ := Ay^^\k = 1,2,3.... If x^^ is sufficient close to 
Xq, then for some constant 7 < 1 

(21) — xoll < — xoll, V/c, 

where 

(22) := argmin{||Q;a:^^^ — Xo|| : |a| = 1, a G C}. 

a 





















Remark 4.2. In view of (16) and Theorem f.l, the following error bound holds for ODR 
with n = N 


where and 

(23) 


- xoll <7^ - Xo||, 


= arg min — a^oll- 

aSC,|o|=l 


Remark 4.3. Theorem 4 A is about the algebraic structure of FDR and does not assume 
oversampled diffraction patterns. For example, one oversampled diffraction pattern (N ^ An) 
or two unoversampled diffraction patterns (N = 2n) are sufficient. 

However, as shown in Theorem 5.3, the proof of (20) reguires one (and only one) over¬ 
sampled coded diffraction pattern. 

Remark 4.4. When the propagation matrix A* is not isometric, we apply QR-decomposition 
to obtain A* = QR, where Q is isometric, and treat Q as the new propagation matrix and 
Rxq as the unknown. 

To this end, we derive and analyze the local approximation of FDR as follows. 

4.1. Local analysis. First note that 

(24) {2l(a<'=V“' - y«). a:o> = («“>!/“' - Vo, 9o> = {V'‘\ Isol) 
with 

- 2 / 0 ), = diag(a;o), 1 x 0 = 

ml 

This motivates the following analysis of the Jacobian operator Jf. 

Proposition 4.5. Let y G u = y/\y\ and D = diag(a;). 

Then 

(25) Sf{y + erf - Si{y) = eklRv + o(e) 
where 

(26) 

with 

(27) 

In particular, if \y\ = b, then (26) becomes 

(28) JfV = {I- B*B)^{v) + iB*B^{v) 

Proof. Let 

w, = ^-r, if = diag(a;J. 

\y + er]\ 

sfy) = y- A*Ay + {2A*A - I)Qb. 


JiV = (J - B*B)v + i{2B*B - I) diag 


B = AVt, v = Ftp. 


[\y\ 


4s{v), 


Reorganizing (18), we have 
(29) 

and hence 


( 30 ) 


Sf{y + ep) - Sf{y) = e{I - A*A)p + {2A*A - I){D, - D)b 

= e(/ - FlB*BFl*)p + {2FIB*BFI* - I){Fl, - Fl)b 













We next give a first order approximation to (f2e — f2)6 in terms of r]. 
Using the first order Taylor expansion we have 


Uf — u = iQQ 


— + o{e) = 


and hence 
(31) 


— r2)6 = ief2diag 


Lbl 


. blJ 

^{n*T]) + oie). 


+ o{e), 


Finally, substituting rtSlh into (30) we obtain 


S{{y + erj) — Si{y) = e(J — VtB*BVt*)r] + ie{2VtB*B — f]) diag(6/|j/|)Q'(r2*r7) + o(e). 
Multiplying fl* on both sides and using the definition of v we complete the proof. □ 

Note that Jf is a real, but not complex, linear map since Ji{cv) ^ cJ^v, c G C in general. 
Dehne the real form of the matrix B\ 


(32) 

Note that 


B : = 


SJ[B1 


e K 


2n,N 


Q[B^] 
^[B^] 

is real isometric because B* is complex isometric. 
From (1^ we have 


(33) 


G{B*u) = 


B^G{u) ' 
B^G{-iu) 


ueC 


Vo 


For the rest of the paper, B denotes the matrix (|27|) with 12 = Qq, i.e. 
(34) 

unless otherwise specihed. 


B = A12o, 12o = diag[a;o], 1:^0 = i—r 

bo I 


4.2. Eigen structure. Let Ai > A2 > ... > A2n > A2n+i = ••• = = 0 be the singular 

values of B with the corresponding right singular vectors {pk ^ and left singular 

vectors {^k ^ By dehnition, for k = 1, ...,2n, 

(35) Br]k = XkG-\^k), 

(36) = ^kVk. 

Proposition 4.6. We have .^1 = G{xo), ^ 2 n = G{—ixo), Ai = 1, A 2 n = 0 as well as rji = |j/o|. 
Proof. Since 

B*xo = ^lA*xo = \yo\ 


we have by (33) 

(37) 

and hence the results. 


3fJ[5*Xo] = = bol, ^B*Xo] =B'f 2 n = 0 


T( 


□ 
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Corollary 4.7. 
(38) 


A2 = max ||m|| ||9=[i?*M]|| = 


u£C 


= max \\u 

ugR^" 
ii_L5i 




Proof. By 

The orthogonality condition iu T Xq is eqnivalent to 

G(xo) T G{—iu). 


Hence, by Proposition 4.6 ^2 is the maximizer of the right hand side of (38), yielding the 
desired value A2. 

□ 


Proposition 4.8. For k = 1,... ,n, 

(39) \\ + \\rL+l-k = 1 

(40) ^2n+i-k = G{—iG ^(^fc)) 

(41) ik = G{iG ^(.^2n+i-fc))- 

Proof. Since B* is an isometry, we have ||tc|| = ||i?*tc||,Vw G C"". On the other hand, we 
have 


and hence 

(42) 


\\B*wf = \\G{B*w)f = \\B^G{w)f + \\B^G{-iw) 


\\GiwW = \\B' G{w)\f + \\B' G{-iw 


ImII = 1 


Now we prove (39), ( |40[ ) and (41) by induction. 

Recall the variational characterization of the singular values/vectors 

(43) Aj = max = argmax s.t. u L 

(44) A2n+l-i = '^i'^ \\^'^u\\, ^2n+l-j = arg miu s.t. U T ^2n, • • . ,6n+2-i, ||ti|| = 1- 

By Proposition 4.6, (39), (40) and (41) hold for k = 1. Suppose (39), (40) and (41) hold for 
/c = 1,..., j — 1 and we now show that they also hold for k = j. 

Hence by (42) 

Aj = max = 1 — min ||H^n||^, s.t. u T ^1 ,..., Q-i, v = G{—iG~^{u)). 


||li||=l 




The condition u T ^1 ,..., P_i implies v T ^ 2 n, ■ ■ ■, ^ 2 n+ 2 -j and vice versa. By (44), we have 


A^ = 1 - Xln+i-j ^^*4 G{-iG H^i)) is the 


mimmizer 


, i.e. ^ 2 n+i-j — G{—iG 


□ 


The relation (39) suggests the following parametrization of singular values of B: 
X 2 n+i-k ■■= COS9k, Afc := sin^fc, 9k G [0,27r]. 


Proposition 4.9. For each k = 1,... ,n, 

(45) B*Brjk = Xk{XkVk + iX 2 n+l-kV 2 n+l-k), 

(46) B Brj2n+l—k •^2n+l—A:(•^2n+l—fch2n+l—fc 
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Proof. By definition, Brjk = Xk^k- Hence 

Btj, = (iR[i?] + t^[B])T], = XkU^ + z^k) 

where 


^k = 


a 


e 


On the other hand, B~^fk = Xkfjk and hence 

(47) Si[B‘"]4'‘ + 

Now we compute B*Br]k as follows. 

(48) = AiB-(4“ + jj;) 

= A^(S|B^a^ + + iA*(»lB^f; - 

= XW. + iA,(Si[sT]5i - 3[ST]4K) 


by 

Notice that 


S(syj'- 3(B' )4» = B 


Cli/'RTAtrR _ 


(49) 


,T 1Xl{-iG-\fk)) 
^-iG-\ik))^ 
= B'^G{-iG-\fk)) 

= i2n+l-k 

•^2n+l—fch2n+l—fe 


by Proposition 4.8 


Putting (48) and (49) together, we have (45). (46) follows from the similar calculation. □ 

Corollary 4.10. For fc = 1, 2,..., 2n, Jf leaves invariant the subspace spanned by {pk, iV 2 n+i-k} 
and has the matrix form 

cos dk sin 9k 


(50) 


J, = A 


2n+l—k 


— sin 9k cos 9k 
in the basis of {pk,iV 2 n+i-k}■ particular, 

(51) Jipi = 0, = ipi 

(52) JiP2n = mn, JpV2n = 0. 
where pi = |?/o|. 


X 2 n+i-k cos9k, Xk ■— sin6'/; 


Proof. By Proposition 4.9, the span of pk and ip 2 n+i-k is invariant under B*B and hence 
under Jf for k = 1, ...,2n. Moreover, (45) and (46) imply 

B*B = 




XkX 


XkX 


k'^2n-\-l—k 


k'^2n+l—k 

X2 

'^2n-\-l-k 


in the basis of i'/72n+i-fc- Hence by the definition (|28|) and Proposition 4.8 

^2n-\-l—k ^k 

^k ^2n+l—k 


Jf = A 


2n+l—fc 


= A 


2n-\-l—k 


cos 9k sin 9k 
— sin 9k cos 9k 


9k € R. 


Hence A2n+i-fe(A2n+i-fc ± iXk) are eigenvalues of Jf. 
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□ 

In the next two propositions, we give a complete characterization of the eigenstructure of 


Ji. 


Proposition 4.11. = 0, A; = 1, 2, 2n — 1, then 

(53) Bv = 0, JfV = 

Proof. The condition v*rik = 0 is equivalent to = rij‘^{v) = 0. So we have 




■ 3?(R5R(n)) ■ 


■ 3fi(R)5R(n) ■ 

S(R5R(n)) 


^{B)ift{v) 


implying B^{v) = 0. Likewise, BQ{v) = 0. So, we have Bv = 
By the dehnition of Jf and Bv = 0, 


= B^{v) = 0 


0 . 


JiV = (J - B*B)^{v) + iB*BQ{v) = 3fi(n). 


□ 


Corollary 4.12. The fixed point set of Jf contains the subspace 

El = nullK(i3) C 

and the null space of Jf contains the subspace 

Eo = iEf. 


Moreover, if X 2 < 1, then 
where 


E 2 — Eq ©r Ef 

E2 = span^{r]kfiVk}l'^f^- 


Proof. Note that 772^ and iri 2 n are excluded from E 2 because r] 2 n G Ei,iri 2 n G Eq. On the 
other hand the null vector pi does not belong in Eq and the eigenvector ipf for eigenvalue 1 

does not belong in Ef for an obvious reason. _ 

For any v G C^, we can write v = 3fJ(n) +i‘^{v). By Proposition 4.11, if 3fJ(n), Q'(n) G Ef, 
then 


Bifi{v) = 0, Jf(3fi(n)) = 5R(n) 

B^{v) = 0, Jf(3(n)) = 0. 


In other words, 3fJ(n) G Ef and 3(n) G Eq. 

On the other hand, if A2 < 1, then A2n-i > 0 and E 2 has no nontrivial intersection with 
either Eq or Ef. Hence, Ef = Eq Ef. 

□ 


4.3. Proof of Theorem 

(54) 


4.1 


With v^^l = -t/o), we can express FDR as 


by Proposition 4.5 Moreover, 

(55) BJfV = {I - BB*)B^{v)+ iBB*B%v) 

= iB^{v), Vn G C^, 


by (28) and the isometry property BB* = L 
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At the optimal phase adjustment for we have 

= 0 

and hence 

(56) {v^^\i\yo\) = - l/oA 2 /o) = * 2 /o) = = 0 

which implies that is orthogonal to the leading right singular vector rji = \yo\ oi B: 

(57) r]l%v^’^'>) =0 


cf. Proposition 4.6 
By (g and 


= ||53(n(^))|| = ||i39(n(^))|| < A2p(n(^))||. 


By induction for fc = 1, 2,... 

(58) ||A(a(''+^) 2 /(^+i) - |/o)il = ||Sn(^+i)|| < 7^||n(^)||, 

for some 7 G [A 2 , 1 ), if is sufficiently small. 

Therefore for = Ay^^\ 

(59) - xo = B - VLly^) = Bv^^\ fc = l,2,... 

which tends to zero geometrically with a rate at most 7 < 1 . 


5 . Spectral gap 


In this section, we prove the spectral gap condition (20) with at least one oversampled 
coded diffraction pattern. This is the immediate consequence of the following two results. 


Proposition 5.1. Let A* be isometric and B = Ahlo- Then ||7s'(i?*x)|| = 1 holds for some 
unit vector x if and only if 

(60) ‘ift{a*x)‘ifi{a*Xo) + 'A{a*x)A{a*Xo) = 0, Vj = 1, ...N, 
where aj are the columns of A, or equivalently 

( 61 ) 717 

where the components of a are either 1 or -1, i.e. 

a(j)e{l,-l}, Vj = l,...iV. 


Proof. We have 


Q{B*x) = 




A*xo 


qA*x 


\A*xo\ 

'iR.{a*xo)'A{a*x) — 7y(a*xo)3f?(a*a:) 
i=i 


(3?2(a*xo) + A2(a*Xo))^/2 
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(62) 









and hence 


N 


||9(i?*a:)|p < + ^'^{a*x) = 

i=i i=i 

by the Canchy-Schwartz ineqnality and isometry. 


N 

Y1 


= ||A*x|p = ||a:| 


In view of (62), the ineqnality becomes an eqnality if and only if (60) or (61) holds. □ 


Proposition 5.2. (Uniqueness of Fourier magnitude retrieval) Let xo be a given rank > 2 
object and /x be continuously and independently distributed on the unit circle. Let 

A*xo 


5 := A diagjtno} , ujq := 


|A*Xo| 


where A* is the given by (§. If 

(63) = ±/3l*a:o 

(after proper adjustment of the angles wherever the coded diffraction patterns vanish) where 
the ± sign may be pixel-dependent, then almost surely x = cxq for some constant c G K. 


The proof of Proposition 5.2 is given in Appendix A. 

Now we can prove the spectral gap theorem needed for geometric convergence of FDR. 

Theorem 5.3. Let $ be the oversampled discrete Fourier transform. Let Xq be a rank > 2 
object and at least one of p.j,j = > 2, be independently and continuously distributed 

on the unit circle. Let 

4 diag{/^i} 


(64) 

be isometric. Let 


A* = c 


$ diag{/X£} 


R := A diagjwo} , Wq : = 




iA*xor 

Then with probability one 

(65) ||A(R*m)|| = 1, ||n|| = 1 iff u = ±ixo/\\xo\\ 

and hence 


( 66 ) 


A 2 = max IIm 


I - 1 II C> 


^{B*u)\\ < 1. 


Proof. Note that the proof of Proposition 5.1 depends only on the fact that A* is isometric 


and hence holds for at least one coded diffraction pattern, oversampled or not. 

Also, the nniqneness theorem. Proposition |5.2[ clearly holds as long as there is at least 
one oversampled coded diffraction pattern. 


Now Proposition 5.1 says that (65) holds if (63) has a nniqne solntion np to a real constant 


and Proposition 5.2 says that (63) indeed has a nniqne solntion np to a real constant. The 
proof is complete. □ 


We have the following corollary from Theorems 4.1 and 5.3 


Corollary 5.4. Let the assumptions of Theorem 5F be satisfied. Then the local, geometric 


convergence (21)-(22) holds for phase retrieval with (64) as the propagation matrix. 
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6 . Uniqueness of fixed point 
Let Doo be a fixed point of FDR ( |I^ , i.e. 

*S'f(|/oo) VoO) ^oo ^Voo- 

Let (Uoo = yoo/lVool be the phase factor of the fixed point. Let 

(67) X = [A{2bQUoc - yoo)]x = [‘^MbQuioo) - Xoo]x ^ 

where X represents the sector condition in the 1-pattern case and df = in the 2-pattern 

A* X = b Q ujoo 


case. We have from (12) 
( 68 ) 

which implies 

(69) 

(70) 


\A*x\ = |A*a:o| 
AA*x = Zi/oo- 


Eq. (69) is related to phase retrieval and eq. (70) magnitude retrieval problem. 
Now we state the uniqueness for Fourier phase retrieval. 


orem 


Proposition 6.1. (Uniqueness of Fourier phase retrieval) Let the assumptions of The- 
5.3 hold. Let x be a solution of of the phase retrieval problem ([^. 


1) One-mask case i = 1. Suppose, in addition, that Axo{j) G [—avr,Vj. Then 
X = e'‘^xo for some constant 0 G K with a high probability which has a simple, lower bound 

(71) 


n 


(3 Sot 


if /i is uniformly distributed on the unit circle, where S is the sparsity of the image and 
[[S'/2JJ the greatest integer less than or equal to S/2. 

2) Two-mask case 7 = 2. Suppose, in addition, that both masks are independently and 
continuously distributed on the unit circle and independent from each other. Then x = e^^xo 
for some constant 6 E R with probability one. 


IS given m 


The proof of Proposition 6.1 
be found, including the l|-mask case. 


27 where more general uniqueness theorems can 


Theorem 6.2. (Uniqueness of fixed point) Under the set-up of Proposition 6.1, the following 
statements hold for the phase retrieval problem ([^. 

1) One-mask case. With probability at least given in (0. X = for some 0 G K. 

2) Two-mask case. Almost surely x = x^o = c^^Xq for some constant 6 E R. 


Remark 6. 3. W ith a slightly stron ger a ssumption ( (211, Theorem 6), the two-mask case in 
Proposition\6.l\ and hence Theorem 6.2 hold for the l^mask case (|4]). 


Proof. By Proposition |6. l| ([6^ implies that x = e^^xo for some constant 9 E [0, 27r], with the 
only difference between case 1) and case 2) being the probability with which this statement 
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Figure 2. The original phantom without phase randomization (left), the 
truncated cameraman (middle) and the truncated Barbara (right). 


holds. To complete the proof, we only need to consider case 2) and show x = Xoo- 

By (67) X G df and by Proposition |6.1| ([^ implies that with probability no less than ( [tI] ) 
X = for some constant 0 G K. Hence, by ([70|, we have 

(72) 


iB 

e 0J{) = 0Jo 


Substituting (72) into (67) we obtain 


2e®^xo = X + Xoo = + Xoo 


and hence e®®Xo = x^o- In other words. 


* 0 , 


Xoo = x = e xo 


□ 


On the other hand, for case 1), x = e®®Xo and (70) imply e®^a;o = oj^o and hence 

QuJo = bQuJoo- 


This and (67) only lead to 


X = [2e®^xo - Xoo]a^ = [2 x - Xoolx- 


7. Numerical simulations 

7.1. Test images. For test images xq we consider the Randomly Phased Phantom (RPP) 
Fig. 1^ (left) and the deterministic image, hereby called the Truncated Cameraman-Barbara 
(TCB), whose real part is the truncated cameraman. Fig. (middle) and whose imaginary 
part is the truncated Barbara, Fig. (right). The purpose of truncation is to create an 
unknown, loose support (dark margins) which makes the image more difficult to recover. 
RPP has a loose support without additional truncation. Likewise, we randomize the original 
phantom in order to make its reconstruction more challenging. In general, a random object 
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Figure 3. Error (on the log-scale) versus iteration for TCB (left) and RPP (right). 


such as RPP is more difficult to recover than a deterministic object such as TCB (see, e.g. 
Fig. I and[^. The size n of both images is 256 x 256, including the margins. 

The propagation matrix is primarily based on either (|^ or (|^ unless specihed otherwise. 

7.2. Local convergence rate. First we simulate the local convergence rate of the 1^-mask 
case and compare them with A 2 . 

The initial condition is chosen sufficiently close to the true object Xq, which is a 
unit vector. Fig. shows the error — Xo|| on the log scale versus the iteration 

counter in the case of two oversampled diffraction patterns. The oscillation in the blue 
curve (FDR) is due to the complex eigenvalues of Jf. The magenta line shows the geometric 
sequence The A 2 value is computed via the power method, A 2 = 0.9505 for TCB 

and A 2 = 0.9533 for RPP. Note that the FDR curve (red) decays slightly faster than the 
A 2 -curve (magenta), which decays still faster than the black curve (ODR with n k. An). 


7.3. Initiailization. For global convergence behaviors, we test two different initializations: 
the Random Initialization (RI), where each pixel value is selected randomly and indepen¬ 
dently, and the Constant Initialization (Cl), where each pixel value is set to unity. 

The relative error of the estimate x with the optimal phase adjustment is given by 


(73) 


ex — Xo 

ll^oll 


je- 


9 = argmin ||e"“x — xnl 
6»eiR 


7.4. One-pattern case. Fig. [^(a)-(d) shows the results of the 1-pattern case for which the 
sector condition is imposed (ODR is equivalent to FDR as N = h). 

To test the effect of the sector constraint, the phase of RPP is uniformly distributed in 
two different intervals: [0,7r/2] and [0,7r]. While FDR/ODR global convergence regardless 
the initialization is evident, the rate of convergence decreases as the sector enlarges. When 
the sector constraint is absent, the iteration ceases to converge in general. 


7.5. l|-mask case. For two- or multi-pattern case, we let the phase of RPP be uniformly 
distributed in [0,27r] (i.e. no sector constraint). Fig. (e)-(h) shows the results of the 1 — 
mask case for which ODR is implemented with n ^ An < N ^ 8n and hence not equivalent 
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(a) RPP + RI 


(b) RPP + Cl 


(c) TCB + RI 


(d) TCB + Cl 



(f) RPP + Cl 


(g) TCB + RI 


— FDR I 

+ odrI 


400 600 800 1000 


(h) TCB + Cl 



Figure 4. Relative error versus iteration in the 1-pattern case with two dif¬ 
ferent sector constraints: [0,7r/2] for (a)-(d) and [0,7r] for (e)-(h). 



Figure 5. Relative error versus iteration in the 1^-mask case. 


to FDR. We see that the performances of ODR and FDR are drastically different: While 
FDR converges to the true images regardless the initialization within 100 iterations, ODR 
does so only for the deterministic image TCB. 

Fig. g shows the relative error versus noise-to-signal ratio (NSR) when noise e is present 
in the data where 


NSR = 


iiki*xoir 

We note that for NSR G [0,20%] there is essentially no difference between the results with 
the maximum number of iterations set to 100 and 200 and this segment of error-noise curves 
is approximately the same straight line (slope ~ 2.2). For a higher NSR, increasing the 
maximum number of iterations reduces the error so with even greater number of iterations 
the straight line segment can be extended to NSR greater than 20%. 
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(b) TCB 


Figure 6 . Relative error versus NSR in the l|-mask case (both oversampled) 
with the maximum number of iteration set to 100 or 200. 
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qI- --- 1 -1 
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(b) RPP + Cl 







Figure 7. Relative error versus iteration with 3 patterns (a)-(d) and 4 pat¬ 
terns (e)-(h) (without oversampling in each pattern). 


7.6. Multi-mask case. To test how DR performs in the setting of multiple patterns with¬ 
out oversampling 14,15 we simulate the 3-pattern and 4-pattern cases with the propagation 
matrices 


(74) 


$ diag{/ii} 


A* =c 


$ diag{/i£_i} ’ 
$ 


£ = 3,4, 


where $ is the standard (unoversampled) discrete Fourier transform. 
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Figure 8 . Relative error for RPP (left) and TCB (right) with various ratios njn. 


Figure shows the result with three patterns (a)-(d) and four patterns (e)-(h), both 
without oversampling, i.e. iV = 3|A^| and iV = 4|A^|, respectively. Note that the number of 
data with four patterns is half of that with 2 oversampled patterns and yet the performance 
of the former is almost as good as that of the latter. 

Going from three patterns (Fig. (a)-(b)) to four patterns (Fig. (e)-(f)) reduces the 
number of iterations by almost an order of magnitude when RPP is the unknown image. 
The case with TCB has less room for improvement. 

7.7. Padding ratio. Finally we test the effect of the padding ratio n/n on the performance 
of ODR. For each n/n G [4, 8], we conduct 50 trials with independent, random initializations 
and average the relative errors. Recall that h/n = 4 is the standard padding rate and at 
n/n = 8 ODR is equivalent to FDR. 

Fig. shows the averaged relative error versus the ratio n/n, demonstrating an effect of 
phase transition which depends on the number of iterations. As the number of iterations 
increases, the threshold ratio decreases. The phase transition accounts for the sharp tran¬ 
sition from stagnation at the standard ratio n/n = 4 to the rapid convergence at n/n = 8 
(i.e. FDR) seen in Fig. (a)-(b) and Fig. 0(a)-(b). (e)-(f) . 

Appendix A. Proof of Proposition 15.21 

In order to prove the uniqueness theorem for Fourier magnitude retrieval, we need to take 
up the more elaborate notation in Section |1.1[ 

Let 

F(z) = 5^/(n)z-" 

n 

be the 2 ;-transform of /. According to the fundamental theorem of algebra, R(z) can be 
written uniquely as 

p 

(75) F{z) = q;z“"o JJ Rfc(z), 

k=l 

where no is a vector of nonnegative integers, a is a complex coefficient, and Ffc(z) are 
nontrivial irreducible polynomials in z~^. 
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Define the shift 


/m+(-) =/(m +■), /m-(-) =/(m-■)• 

Definition 1 (Conjugate Symmetry). A polynomial X{z) in z~^ is said to be conjugate 
symmetric if, for some vector k of positive integers and some 9 G [0,27r), 

X{z) = 

In other words, the ratio between X{z) and its conjugate inversion is a monomial in z~^ 
times a complex number of unit modulus. 

A conjugate symmetric polynomial may be reducible, irreducible, trivial, or nontrivial. 
Any monomial z~^ is conjugate symmetric. 

Proposition A.l. Let f be a finite array whose z-transform has no conjugate symmetric 
factors. If the z-transform G of another array g satisfies G {0,7r}, 

Vw G C then g = cf for some constant c G K. 

The real-valued version of the above proposition is given in [^. For the reader’s conve¬ 
nience, we provide the proof for the complex setting. 


Proof. Consider the array h dehned by 

h = f-k g{—) 

whose 2 ;-transform is 

^ H{z) = F(z)G(F^. 

Note that h is dehned on AA, instead of Ad, so hf(z) is completely determined by sampling 
H on C. 

Since 

it follows is real-valued. By analytic continuation, we have 

H{z) = id(z-i) 

and 


(76) 


F(z)G(z-i) = F(z-i)G(z). 


Multiplying both sides of (76) by z ^ results in the following polynomial equation in z 


(77) 


F(z)G'(z 


-n 


-M 


= Z 


-M 


F(z-i)G(z). 


We observe no = 0 in view of (75) and the assumption that F{z) has no conjugate 
symmetric factor. We also have 

(78) z~'^F(z“^) = az~ 

k 




where Fk{z) are the nontrivial irreducible non-conjugate symmetric polynomials in z ^ of 
the form Fk{z) = Fk{z~^) for some vector of positive integers. 

Writing 

(79) G(z) = 0z-”»nGf(z), 

I 
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where Ge{z) are nontrivial irreducible polynomials in z we have 
(80) = /3z—1 H G,(z), 


where Gi{z) are the nontrivial irreducible polynomials in z~^ of the form G£{z) = z~^~^'^^Gf (z~^^ 
for some vector of positive integers. 

and 


Plugging (|75|),(|78|), 
(81) 


in (77) yields 


a/3z—1 H F,{z) J] Ge{z) = d/3z-—« J] Ffc(z) J] Ge{z). 


Each nontrivial irreducible factor Efc(z) must be equal to some Fk/{z) or some Gii{z). How¬ 
ever, if iq.(z) = Fk{z), then ihc(z) is a conjugate symmetric factor. If, on the other hand, 
Ffc(z) = Fk'iz) for some k’ ^ k, then Fk{z)Fk'{z) = Fki{z)Fk{z) is a conjugate symmetric 
factor. Both cases, however, are excluded by the assumption that the ; 2 -transform of / does 
not have conjugate symmetric factors. 

Hence each Fk (rest. Fk) must be equal to some G^ (rest. Gg) and we can write 
(82) G'(z) = Q(z)F(z) 

where Q{z) is a polynomial in z~^, i.e. 

Q{z) = y^CnZ~". 

n>0 

By the assumption that ZF(e^”'^) — ZG'(e^’^*’^) G {0, vr} we have (^(e^”"^) G [R,Vw G £, 
and hence Cn = c_n = 0 except for n = 0 in which case cq G K. Therefore, Q = cq G K and 
this is what we start out to prove. □ 


Proposition A.2. Let xq have rank > 2. Let {/i(n)} he independent and eontinuous 
random variables on the unit circle of the complex plane. Then, the z-transform F{z) of 
/(n) := /i(n)xo(n) is irreducible up to a power of z~^ with probability one. 


For the proof of Proposition A.2 see Theorem 2 of 27 . 


We next show that the 2 ;-transform of {p(n)a;o(n)} is almost surely irreducible up to a 
power z~^ and not conjugate symmetric. 


Proposition A.3. Let {p(n)} he independent and continuous random variables on the unit 
circle of the complex plane. Let /(n) := /r(n)a;o(n). Then the z-transforms of both /t+ and 
ft- are almost surely not conjugate symmetric V t. 


^t+(z) = ^/(t + n)z 


Proof. The ^-transform 

(83) 

is conjugate symmetric if 

(84) 

for some vector k of positive integers and some 9 G [0,27r). Plugging (83) in (84) yields 


Ft+(z) = e*^z-*^Ft+(z-i) 


fit + n)z-’^ = F^z-^ 
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which implies 

(85) /(t + n) = e*®/(t + k — n), Vn. 


However, xo is deterministic, and {p(n)} are independent and continnous random variables 
on S^, so (85) fails with probability one for any k. There are finitely many choices of k, so 
the ^-transform of /t+ is almost snrely not conjngate symmetric. 

Similarly, the ^-transform of ft- is also almost snrely not conjngate symmetric. □ 
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